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Abstract 

In this paper we show that the total variation diminishing (TVD) finite 
difference scheme which was analysed by Sweby [6] can be interpreted as a Lax- 
Wendroff scheme plus an upwind weighted artificial dissipation term. We then 
show that if we choose a particular flux limiter and remove the requirement 
for upwind weighting, we obtain an artificial dissipation term which is based 
on the theory of TVD schemes, which does not contain any problem dependent 
parameters and which can be added to existing MacCormack method codes. 
Finally, we conduct numerical experiments to examine the performance of this 
new method. 
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1* Introduction 


A major result of recent research into numerical methods for the solution 
of systems of hyperbolic conservation laws has been the development of second 
order accurate total variation diminishing (TVD) finite difference schemes. 
These schemes have a number of attractive properties including the fact that 
they resolve discontinuous solutions well, they do not exhibit spurious 
oscillations and, under certain circumstances, they can be proven to converge. 

Unfortunately, the apparent complexity of these schemes has thus far 
discouraged their widespread adoption. Instead, most applications progams use 
standard methods such as the MacCormack variant of the Lax-Wendroff method and 
add additional artificial dissipation to damp the spurious wiggles that occur 
near discontinuities. These dissipation terms are usually chosen in an ad hoc 
fashion and usually contain problem-dependent parameters which must be fine 
tuned before the method will work. 

In this paper we attempt to remedy this situation. In particular, we 
interpret the TVD finite difference scheme which was analysed by Sweby [6] as 
a Lax-Wendroff scheme plus an upwind weighted artificial dissipation term. We 
then attempt to simplify this artificial dissipation term by removing the 
requirement that it be upwind weighted. The result is a dissipation term which 
is based on the theory of TVD schemes and which does not contain any free 
parameters. 

An outline of this paper is as follows. In section 2 we briefly review 
the theory of TVD finite difference schemes. Using this theory, we derive an 
artificial dissipation term for scalar hyperbolic equations in section 3 and 
we examine its performance in numerical experiments. These results are 
extended to systems of hyperbolic equations in section 4 and additional 
numerical experiments are performed. Section 5 contains some closing remarks. 
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2. Total Variation Diminishing Finite Difference Schemes 

We consider the initial value problem for a scalar conservation law. That 
is 

u t + f(u) x = u fc - a(u)u x = 0, a(u) = (u), t > 0 

( 2 . 1 ) 

u(x,0) = Uq(x) , -« < x < 09 

where u Q (x) is assumed to have bounded total variation. A weak solution to 
this problem has the following monotonicity properties. 

(1) No new extrema in x may be created. 

(2) The value of a local minimum is nondecreasing and the value of a local 
maximum in nonincreasing. 

The total variation of the solution to (2.1) at time t is defined by the 
formula 

(2.2) TV( u(x, t) ) - sup l |u(x k+1 ,t) - u(x k ',t) | , 

k 

where the supremum is taken over all partitions of the real line. 

It follows from this monotonicity property that the total variation in 
x of u(x,t) does not increase in t. That is 

(2.3) TV(u(t 2 )) _< TV ( u (t 1 )), for all t £ > tj. 

Much recent research has been devoted to the construction of finite 
difference schemes that satisfy a discrete version of equation (2.2). We 
briefly describe this work below. 



Consider explicit finite difference schemes in conservation form which 


approximate (2.1) and which we denote by 
(2.4) U n+1 - L*U n . 


A scheme is called total variation diminishing if 
(2.5) TV(u n+1 ) - TV(L*U n ) < TV(U n ). 


In addition, a scheme is called monotonicity preserving if the finite 
difference operator L is monotonicity preserving; that is, U n a monotone 
mesh function implies that L»U n is also a monotone mesh function. The 
following result, presented without proof is due to Harten [3]. 

Theorem 2.1 (Harten). A total variation diminishing scheme is 
monotonicity preserving. 

This says that TVD schemes will not produce spurious oscillations. This 
is their chief attraction. 

Another reason why TVD schemes are attractive is that it is very useful to 
have a bound on the total variation of the solution when proving convergence 
of nonlinear difference schemes [cf. 4, 5]. Equation (2.5) provides such a 
bound but convergence proofs are beyond the scope of this paper. 

The scheme (2.4) can be written in the form 


( 2 . 6 ) 


"f 1 - d £ - c k-v 2 4 D £-v 2 + WC1/2 
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where 

< 2 - 7 > < + v 2 ‘ ' u " 

and C k _ 1 /2 and D k+1 / 2 are functions of U n . Harten [3] proves the 
following Lemma which provides a sufficient condition for the scheme (2*6) to 
be total variation diminishing. 

Lemma 2.2 (Harten). If the coefficients C and D of equation (2.6) 
satisfy the inequalities 

< 2 - 8 > 0 — D k+ V2 

0<C k+1 ^ + D k+l4 <! 

then the scheme (2.6) is total variation diminishing. 

Later, we use this Lemma to prove that our scheme is TVD. 


3. Scalar Equations 

For ease of presentation, we first consider the scalar linear equation 


(3.1) 


u t + au x 


0, a = const. > 0 
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Sweby [6] considers the solution to this problem using a scheme of the 


form 


(3.2) u£ +1 = u£ - v{l +V 2 (l-v)[<f»(rp/rJ - r+J] }Au£_ ^ 


+\ # + 


where 


(3.3) 


,-ass. 

At k ATT n 

v 2 


and < t , ( r j c ) is the flux limiter. This is a scheme of the form (2.6) with 


(3.4) 


C k _ l /2 - v{l + l / 2 (l-v)[*(r£)/r£ - 


D. 


k+ V 2 ~ 0 * 


A sufficient condition for the scheme (3.2) to be TVD is that v < 1 and 


(3.5) 


!■♦■( r k ) / r k - I < 2. 


Sweby also specifies that <|)(r) > 0 and that 4>(r) = 0 for r 0. Under 
these additional restrictions the bound (3.5) becomes 


(3.6) 


0 < <|)(r)/r, <t> (r) < 2. 


If <|>(r) * 1, scheme (3.2) 
Wendroff method. If ( r) * r, 
upwind Warming and Beam [8] method. 


reduces to the centered difference Lax- 
scheme (3.2) reduces to the second order 
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The region defined by (3.6) is shown in Figure 1 along with the limiters 
corresponding to the Lax-Wendroff and Warming-Beam methods. Since these 
schemes are known to produce spurious wiggles in solutions with steep 
gradients, it is not surprising that these schemes are not uniformly within 
the TVD region. 

Since the Lax-Wendroff method does not require one to determine an upwind 
direction and since many production computer codes are based on the Lax- 
Wendroff method or its variants, we wish to examine the possibility of adding 
terms to these codes to obtain a TVD scheme. 

If we add a term of the form 


(3.7) 


v 2 ( ■><£. 1 /2 - <_ i /2 ( ^K- 


to the Lax-Wendroff scheme 


(3.8) uf 1 - B” l/ 2 ) + f K+V 4B k-V 2 ) 


and rearrange this to the form (3.2) we get 


(3.9) Uj +1 = Uj - [v{l +V 2 (1 -v)[1/4 - 1]} " K^_ 1/2 }]AU^_ 1/2 


That is 


<Vi /2 - v!l + V 2 d-v)[i/4 - i]} - 


D. . i, = 0. 


k+ V 2 


(3.10) 
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A comparison of (3.9) with (3.2) shows that we can obtain Sweby's scheme 
if we choose 

(3.ii) 4+vi a_v) i i - +(>$]• 

Next we consider the equation 


(3.12) 


u t + au x = 0, a = const. < 0. 


For this problem Sweby's scheme takes the form 


♦( r J • 


(3.13) U£ +1 = u£ + v{-l + V 2 (l+v)[$(r~ +1 ) ^-]} AU k + V 2 


where 


au: 


n 


r k = 


k+1 /2 m = jAt 


AU 


n 


Ax 


< 0 


k- V 2 


and our modified Lax-Wendroff scheme takes the form 


(3.14) u £ +1 - U£ + v[{-l + V 2 (l+v)[ 1 } + {\ + ! /2 -%^}]au“ + i /2 . 


L k 


This reduces to Sweby's scheme if we choose 


(3.15) 




A comparison of equations (3.11) and (3.15) shows how the sign of the 
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coefficient a changes the method. In particular, we see that as a changes 
sign, the definition of r changes and the form of the term involving the 
Courant number v changes. We can combine these two cases into a Lax- 
Wendroff method with an upstream weighted artificial dissipation term as 
follows. Put the scheme into the form 


(3.16) 



\ n - 1 tCi- Ci) + 4 ti, - K + 

+ [>W r D + 

- t<-v 2 ( r k-ii + VlyS,( r k)]( U k-°k-l) 


where 


(3.17) 



J (l-v)[l - *(<)] 
0 


if a > 0 


if’ a < 0 



0 

j (l+v)|>(r" +1 ) - 1] 


if a > 0 
if a < 0. 


This method still requires that we know which direction is upwind. For 
hyperbolic systems it is this requirement that makes upwind difference schemes 
complicated. Therefore we attempt to construct a method without upstream 
weighting by rewriting equation (3.17) as follows 
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J £ + V J T L(1 -M>[1 -*(«$] 

(3.18) 

^+1/2“ (1 "M>U “^ r k+ l^* 

It is obvious that the terms involving the K's are dissipative. The 
question that needs to be adressed is how much more dissipative (3.18) is than 
(3.17). To obtain some idea about this we compute the solution to the 
equation 

(3.19) u t + u x = 0 

with square wave initial data. We use the MacCormack version of the Lax- 
Wendroff method and the limiter defined by 

( min(2r, 1) , if r > 0 

(3.20) <f>(r) = 

(0 , if r < 0. 

Figure 2a shows the result of this computation after 100 steps at a 
Courant number, v = .9 using the MacCormack method without additional 
dissipation. Note that the solution exhibits severe oscillations in those 
regions where r + is small as predicted by Figure 1. Figure 2b is the result 
of the same computation using the MacCormack scheme with the upwind 
dissipation (3.17) and the flux limiter (3.20). Notice that the spurious 
oscillations have been removed. Finally, Figure 2c shows the result of the 
same computation using the MacCormack scheme with the simplified dissipation 
(3.18) and the flux limiter (3.20). These results are almost indistinguish- 
able from those of Figure 2b. In addition we can prove the following: 



- 10 - 


Theorem 3.1. The method (3.16), (3.18) with flux limiter (3.20) is TVD 
under the restriction that the Courant number |v| < 1. 

Proof. The proof is a direct appliction of Lemma 2.2. That is, we show 
that if the scheme is put into the form 

(3.21) U " +1 - Ujj - C k _ 1/24 u£_ 1/2+ \ +1/2 M" +V2 

then 


(3.22) 


"k+ V2 


> 0 


(3.23) 


\+V 2 


(3.24) 


0 < V i /2 + \ + i k < i- 


Here we prove the result for v > 0. The computations for v < 0 are 
similar. 

Rearrange (3.16) as follows 


(3.25) 


of 1 Mi+Vza-vXV 1 )! 

r k 

- ( Kk + 1 <2_ „+ 


(-M- <-i 4 - Vfe'K-i** w b ^v 2 


That is 
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<3-26) C fc+ 1/ * v{l +1* <l-v)(iy -1)) - f^2Zi - 4 + 1/2 - 1/2 1 


k+1 


k+1 


(3.27) 


Dfc+ V 2 *k+ V 2 * 


\ + y 2 >° V%> °- 

To show that C ^ + y > 0, substitute (3.18) into (3.26) and rearrange terms 


to obtain 


(3.28) 


: K+ l/ 2 


. v(l-v) f^k+l) , , 

= v + +(r ) +1 - 


^k+l^ 


k+1 


Note that (3.20) implies 


(3.29) 


0 < Mil < 2 

r 


(3.30) 


0 < <f>(r) < 1 


so 


<W/ 2 >v-^-* + 4>0. 


Since %+ l/ 2 > 0 and C k+ 1 ^ > 0, C k+ l/ 2 + D k+ 1 0. To prove that 


C k+ V 2 + °k+ V 2 < 1 ’ 


substitute (3.18) into (3.27) and add to (3.28). The result is 


C3-31) V 1/2 + D fc+ l/2 - V + - ♦( r*) + 2 - 2*( r^,)] 


k+1 
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(3.31) is < 1, if the term in the square brackets is < 2. To show this, 
we note that 


(3.32) 

and we consider four cases. 


r k+l + 
r k+l 


5®as-L- 4+1 < 0 => -r- < 0 => K4+i) = 

r k+l 

In this case the term in square brackets becomes 


- 0 . 

r k+l 


2 - ■♦( 4 ) < 2 


Case 2. 2r£ +1 < 1 => rj +1 < l / 2 > >2, so 

r k+l 

^k+l) = 2r k+l ’ ‘K-*—) " 1 

r k+l 

and the term in square brackets is 


2 - KrJ < 2. 


Case 3. V 2 < Vl' ^- < 2 -> - 1. 


k+1 


k+1 


Then the term in square brackets becomes 


-J_ - <0(4) < 2 - <0(4) < 2. 

r k+l 
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Case 4. 


2 < r* i 

— k+l 


s> ~r~ < 1 => ^ r k+i^ “ l > ^-r~) =-r 


k+l 


“k+l k+l 


Then the term in square brackets becomes 


■+h - *(4) + 2 - --±- - 2 - - *(r£ < 2. 


k+l 


‘k+l 


k+l 


This completes the proof. 


It is trivial to extend these results to scalar nonlinear problems. We 
simply define a local wave speed by 


(3.33) 


\+l/2 



i£ 4 \ + i / 2 


* 0 


11 4D k + V2-0 


and apply the schemes (3.16), (3.17) or (3.16), (3.18) as before. The wave 
speed definition (3.33) makes the resulting scheme conservative. 

Figures 3a, 3b and 3c show computed solutions of the inviscid Burgers' 
equation with square wave initial data and periodic boundary conditions. 
These results were obtained using the MacCormack scheme, the second order 
upwind scheme (3.16), (3.17) and the simplified scheme (3.16), (3.18), 

respectively. The MacCormack results exhibit severe oscillations in the 
vicinity of the shock and an entropy violating expansion shock. The upwind 
scheme eliminates the oscillations but not the expansion shock while the 
simplified scheme eliminates all but a small "entropy glitch" in the expansion 
region. There is no difference between the upwind method and the simplified 
method in their ability to resolve the shock. 
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The Burgers' equation results were highly dependent on the choice of 
initial conditions. Indeed, for certain initial conditions, the simplified 
method also conoputed solutions containing expansion shocks. To our knowledge, 
the only way to avoid expansion shocks in all cases is to explicitly add extra 
dissipation to the method when a sonic point occurs in an expansion region. 
Although we cannot guarantee that the simplified scheme will do this 
automatically, our computations indicate that this scheme is more robust than 
the unmodified upwind scheme. 


4. Hyperbolic Systems 

In this section we extend the results of the previous section to 
hyperbolic systems. To that end, we consider first the linear, constant 
coefficient system 

(4.1) u fc + Au x “ 0, A = const. 

where u is an m vector and A is an m x m matrix. 

If the system (4.1) is hyperbolic, the matrix A has real eigenvalues and 
a complete set of linearly independent right eigenvectors. If we let P 
denote the matrix whose columns are the right eigenvectors of A, then 

(4.2) P -1 AP = A 


where 
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(4.3) 


and are the eigenvalues of A. 


m -J 


If we define a new set of dependent variables by the formula 


(4.4) 


D -1 

v = P u 


and multiply (4.1) by P”* we obtain 

(4.5) (P^u) + P -1 AP _1 (u) = 0 

t X 

or 

(4.6) v. + Av = 0. 

N ' t x 

This is an uncoupled set of scalar equations. We solve (4.6) using (3.16). 
That is 


(4.7) 



where v = AAt/Ax and K* and r* are defined below. 

Multiply (4.7) by P to obtain an equation in terms of the original 
dependent variables. The result is 
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n+1 

J k 


U, n - A 


2Ax 


C<« - <-i) J 


(4.8) 


+ p t< + i/ 2 «) + V^'miD'X. " D k) 


The first three terms on the right of (4.8) comprise the well-known Lax- 
Wendroff scheme and need no further discussion. Hie last two terms require 
that we know the matrices P and P”^ and that we know which direction is 
upwind. In the following we construct a simplified version of this scheme 
which removes these requirements. 

We remove the requirement that P and P"”* be known by approximating the 
diagonal matrices K* by scalar matrices. That is, we let 


(4.9) 


K^r*) « K ± (r ± ) I 


where the K [r J are scalar functions of r . 

We remove the necessity to determine upwind directions by choosing 


(4.10) 


K *(r*) = .5 C(v) [l - Kr*)] 


where the Courant number v is defined as 


(4.11) 


l, | At 
v = ma x X . t— 

3 J 4x 


and C(v) is chosen as follows 
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(4.12) 


C(v) 


v(l-v), 
, .25, 


if v < .5 
if v > .5 , 


This definition of C(v) is an upper bound to the Courant number- 
dependent coefficient in (3.18). 

Thus far we have not defined r + and r~. In light of the fact that we 
do not wish to conpute P or P“^, we have chosen the following definitions 
for r + and r“. 


(4. 1 3a) 


(4.13b) 


k K + v<v 2 ) 
k K. v «Cv 2 ) 


where (•»•) denotes the usual inner product on R m . 

If P does not vary significantly over adjacent mesh intervals, these 
definitions can be interpreted as averages of the scalar definitions. These 
were the most simple definitions that we could construct. They worked so well 
in our numerical experiments that we saw no reason to investigate more 
sophisticated r* definitions. We note in passing that other definitions 
of r* have been proposed by Sweby [6] and Chakravarthy and Os her [2]. 

With these simplifications, our numerical method takes the form 
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TT n+l TT n At f IT n n >> .1 At , n n n ^ 

\ ■ °k ‘ A in (Vi ’ ViJ + A 2 ET l D k-H ' 2 \ + Vi) 


.2 At , n 


(4.14) 


+ + \; 1/2 t r ; +1 )][uj +1 - -op 


- (Cv 2 (vi) + v'lii'PiK - Ci) 


where K and r are defined by equations (4.10) and (4.13) 
respectively. Note that the resulting scheme does not depend explicitly on 
the transformation (4.4). Therefore, we can use the scheme without 
modification on nonlinear problems where (4.5) is not true. For the 
computations which follow, we replace the Lax-Wendroff scheme (the first three 
terms on the right of (4.14)) by the conservative MacCormack scheme. These 
schemes are equivalent for linear problems. 

As a first test, we demonstrate the performance of this method on the 
Riemann problem. That is, we solve the Euler equations 


(4.15a) 


u + f(u) =0, -“ < x < “, t > 0 

t x 


where 




" p " 


- pu 

(4.15b) 

u - 

m 

- E - 

, f(u) = 

2 . 

pu + p 
. (E+p)u _ 


(4.15c) 


p = (y-1)(E - V 2 Pu 2 ) 


and Y - 1,4 with initial conditions 
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(4.15d) 


u(x,0) 



x < 0 
x > 0 


In this case the initial conditions were 


(4.16) 



Figure 4 shows the solution after 100 time steps computed on 140 grid 
points at a Courant number of .95. These conditions are the same as those 
used by Harten [3] and Chakravarty and Qsher [2] in their numerical 
experiments. The results shown in Figure 4 cannot be distinguished from those 
shown in the cited references with one exception. Harten was able to obtain a 
dramatic improvement in the resolution of the contact discontinuity by 
selectively adding artificial compression to his second order upwind scheme. 
We intend to study this technique in the future. 

Next we demonstrate the performance of our method on a two dimensional 
problem. Figure 5 shows a comparison of the present method with the second 
order upwind scheme of van Leer [7], To obtain these results, we solve the 
two dimensional Euler equations 


(4.17a) 


u. + f(u) + g(u) '« 0 
t ' x e y 


where 


(4.17b) 


- p - 


r p u " 


pv 



2 , 



pu 


pu + p 


puv 

pv 

, f(u) = 


, g(u) * 

2 

puv 

pv + p 

- E . 


_ (E+p)u - 


- (E+p)v- 
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(4,17c) p = (y-1) [e - V 2 (pu 2 + pv 2 )] 

and Y = 1.4 for the problem of the reflection of an oblique shock from a 
plane wall. 

For the computations shown, we specify a uniform M ** 2.9 flow at the 
left boundary. At the top boundary, we specify the conditions behind a shock 
that would turn the flow 11° (cf. [1]). A flow tangency condition is 
specified at the wall and all variables are extrapolated at the right 
boundary. The computation is started with the upstream conditions specified 
everywhere except the top boundary. The results shown consist of a three 
dimensional plot of the converged density solution and a longitudinal section 
of this plot taken at y - .5. Once again it is difficult to distinguish 
between the results computed using the two methods. 

Finally, we compare the present method to the second order upwind method 
of van Leer on a model transonic flow problem. 

The two dimensional Euler equations (4.17) are solved for the flow over a 
10% thick parabolic arc bump in a channel. The flow is assumed to be uniform 
initially and then the condition that the flow be tangent to the bump is 
applied at the wall in the manner of small disturbance theory. Nonreflecting 
boundary conditions are applied at both upstream and downstream boundaries. 
Figures 6a and 6b show the converged pressure distribution on the wall for a 
flow with inlet Mach number .675 using the van Leer method and the present 
method, respectively. 

The flow consists of an expansion region over the forward part of the 
bump, 0 £ x < .5, followed by a shock near x = .75. This is the type of 
problem that the Van Leer method was designed to solve since the shock is 
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steady and nearly aligned with the computing grid. Figure 6a shows that the 
van Leer scheme does an admirable job on this problem. In particular, the 
flow expands smoothly through the sonic point (P « .717) and there is only 
one grid point within the shock. By comparison, Figure 6b shows that the flow 
computed using the present method is not quite as smooth in the expansion 
region and there are two points within the shock. Still, considering the fact 
that the current method is considerably easier to program and that it runs in 
2/3 the time of the Van Leer scheme, the results shown are quite acceptable. 

During the course of this work, we discovered that, for two dimensional 
problems, the flux limiter (3.20) imposes a severe Courant number stability 
restriction on the method beyond that of the two-dimensional MacCormack 
scheme. To prevent that we define a new flux limiter by the formula 

(4.18) <p ( r ) - min(2|r| , 1), 

A simple application of the maximum principle shows that the artificial 
dissipation based on this limiter is stable for two dimensional problems under 
the same Courant number restriction as the two-dimensional MacCormack 
method. At this time we have not analysed the three dimensional case. 


5. Concluding Remarks 

In this paper we construct a simple artificial viscosity term for Lax- 
Wendroff type methods which is based on the theory of Total Variation 
Diminishing upwind finite difference schemes. This method has advantages over 
both conventional artificial viscosity schemes and the TVD upwind schemes. In 
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particular, the method does not contain the problem dependent parameters of 
conventional artificial viscosity schemes and it does not require the complex 
logic of upwind schemes. 

The numerical experiments that have been performed thus far have been very 
encouraging but more numerical experimentation is needed. We intend to carry 
this out in the future and also to apply the ideas of this paper to other 
numerical schemes. 
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Figure 1. TVD Region 




Figure 2a. Solution of (3.19) with square wave initial data after 
100 steps at v ■» .9 using MacCormack method. 




Figure 2b. Solution of (3.19) with square wave initial data after 

100 steps at v = .9 using upwind scheme (3.16), (3.17). 




Figure 2c. Solution of (3.19) with square wave initial data after 

100 time steps at v = .9 using TVD scheme (3.16), (3.18). 
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Solution to Riemann problem after 100 
v = .95 using scheme (4.14), (4.10) . 


time steps a 
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Figure 6b. Pressure distribution along wall for transonic problem using 


modified MacCormack scheme. 
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